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ABSTRACT 

We report results from numerical simulations of star formation in the early universe that focus on 
gas at very high densities and very low metallicities. We argue that the gas in the central regions of 
protogalactic halos will fragment as long as it carries sufficient angular momentum. Rotation leads 
to the build-up of massive disk-like structures which fragment to form protostars. At metallicities 
Z k 10~ 5 Zq, dust cooling becomes effective and leads to a sudden drop of temperature at densities 
above n — 10 12 cm~ 3 . This induces vigorous fragmentation, leading to a very densely-packed cluster 
of low- mass stars. This is the first stellar cluster. The mass function of stars peaks below 1M Q , 
similar to what is found in the solar neighborhood, and comparable to the masses of the very-low 
metallicity subgiant stars recently discovered in the halo of our Milky Way. We find that even purely 
primordial gas can fragment at densities 10 14 cm~ 3 < n < 10 16 cm -3 , although the resulting mass 
function contains only a few objects (at least a factor of ten less than the Z = 1CP 5 Zq mass function), 
and is biased towards higher masses. A similar result is found for gas with Z = 10 Z Q . Gas with Z 
< 10~ 6 Zq behaves roughly isothermally at these densities (with polytropic exponent 7 w 1.06) and 
the massive disk-like structures that form due to angular momentum conservation will be marginally 
unstable. As fragmentation is less efficient, we expect stars with Z < 10 -6 Z Q to be massive, with 
masses in excess of several tens of solar masses, consistent with the results from previous studies. 
Subject headings: stars: formation - stars: mass function - early universe - hydrodynamics - equation 
of state - methods: numerical 



1. INTRODUCTION 

Detailed numerical simulations of the formation of 
the first stars, the so-called population III stars, in- 
dicate that they were probably massive, with masses 



carbo n and oxygen b e comes effective (Bromm et al.l 
200 lj E romm fe Loebl [20031: ISantoro fc Shulll [20061 



aicate tnat tney were propapiy massive, witn masses 
greater than 20 Mq (lAbel. Bryan, fc Normanl 1 2002 



Bromm. Coppi, fc Larsonl 120021 : lYoshida et ali 2006 
O'Shea fc Normanl 12007m ' The fact that no popula- 



tion III star has ever been observed in the Milky Way 
provides some ob servational support for this prediction 
(jTumlinsonl 120061 ). However, a number of low-mass, 
extremely metal-poor, stars with [ Ee/H] < —3 have 
been discovered in the Galactic halo (|Beers fc Christliebl 
|2005| ) , suggesting that the distribution of stellar masses is 
sensitive to even very low levels of metal enrichment. Ex- 
planations of this apparent change in the IMF have con- 
centrated on the fact that metal-enriched gas has more 
coolants than its primordial counterpart. These coolants 
are suggested to provide an opportunity for efficient frag- 
mentation, since they can keep the gas temperature lower 
during the collapse process than is possible with pure H2 
cooling in the primordial gas. If this is the case, then 
the final IMF will likely contain at least some low-mass 
stars, even if it still differs significantly from the present- 
day local IMF. 

If metal enrichment is the key to the formation of 
low-mass stars, then logically there must be some crit- 
ical metallicity Z cr ;t at which the formation of low 
mass stars first becomes possible. However, the value 
of Z cr it is a matter of ongoing debate. Some mod- 
els suggest that low mass star formation becomes pos- 
sible only once atomic fine-structure line cooling from 



Frebel. Johnson, fc Bromml 120071 ). setting a value for 
Z cr it at around 10~ 35 Zq. Another possibility, and the 
one that we explore with this paper, is that low mass 
star formation is a result of dust-induced fragmenta- 
tion occurring at high densities, and thus at a very 
late s tage in the protostellar collapse jS chneidc r~et al.l 
l27MlQ^ru"ikai et alJl2005t ISchneider etal\\2(M) . In this 
model, 10~ 6 < Z crit < 10~ 4 Z Q , where much of the un- 
certainty in the predicted value results from uncertain- 
ties in the dust composition and the d egree of gas-phase 
depletion (jSchneider et al.ll2002l [20061 ). 

The recent simu lations performed by 

iTsuribe fc Omukail ([20061 ) . which model the col- 
lapse of very high density protogalactic gas, provide 
some support for the dust-induced fragmentation model. 
Using a simple piecewise polytropic equation of state to 
describe the ther mal evolution of extremely metal-poor 
protogalactic gas, ITsuribe fc Omukail (|2006h show that 
fragmentation can occur at metallicities as low as 
Z = 10 -6 Zq, and that it becomes more effective as the 
metallicity increases. However, their study considered 
only the limiting case of gas with zero angular momen- 
tum. Owing to the absence of angular momentum, the 
fragments formed in their simulation do not survive for 
longer than a dynamical time, as they simply fall to the 
center of the potential well, where they merge with other 
fragments. It is also unclear whether the fragmentation 
process would be as effective if the angular momentum 
of the gas were non-zero, as would be expected in reality. 
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Fig. 1. — The three equations of state (EOSs) from Omukai et al. 
(2005) that are used in our study. The primordial case (solid line), 
Z = 10 -6 Z Q (dotted line), and Z = 10 -5 Z Q (dashed line), are 
shown alongside an example of a polytropic EOS with an effective 
7 = 1.06. 



We present the results of simulations of the high- 
density, dust-cooling dominated re gime that improve on 
those of iTsuribe fc Omukail (|2006f ) by including the ef- 
fects of rotation, and by following a much larger dy- 
namical range of the collapse (ten orders of magnitude 
in density), as well as employing an equation of state 
(EOS) which follows that of Omukai et al. (2005) more 
closely. We also perform some simulations of purely pri- 
mordial gas for comparison. A key feature is the use 
of sink particles (Bate, Bonnell & Price 1995) to cap- 
ture the formation and evolution of multiple collapsing 
cores, which enables us to follow the evolution of the star- 
forming gas over several free-fall timescales and thus to 
model the build-up of a stellar cluster. This differs from 
previous calculations which eit her follow the collapse of 
a single core to high densit i es (lAbel. Bryan, fc Normanl 
I2002t iBromm fc Loebll200l lYoshida et al\\200$) , or use 



sink particles to c apture low density (n < 10 6 c m 3 ) 
fragmentation (e.g. IBromm. Coppi. fc Larsonll2002f ). 

In following section, we give details of the numerical 
model (<J2J} and the initial set-up of the simulations 
( §2.2p . The main results of our study are presented in 
U21 and we discuss the origin of the fragmentation in SJH 
Potential caveats with the current model are highlighted 
in 53 and there is a summary of our findings in ||6] 

2. DETAILS OF THE CALCULATIONS 
2.1. Numerical Method 

We follow the evolution of the gas in this study us- 
ing Smoothed Particle Hydrodynamics (SPH). Our ver- 
sion of the code is essentially a parallelized version of 
that used by Bate et al. (1995), which utilizes adap- 
tive smoothing lengths for the gas particles and a bi- 
nary tree algorithm for computing gravitational forces 
and constructing SPH neighbor lists. The calculations 
were performed on the Julich Multi-processor (JUMP) 
supercomputer at the John von Neumann Institute for 
Computing, Research Center Julich, Germany 

The thermal evolution of the gas in our simulations is 
modelled using a tabulated equation of state (EOS) that 
is based on the r e sults of the detailed chemical modelling 
of lOmukai etlll (12005). We use their reported results for 
the temperature and molecular fraction as functions of 
gas density (their Figures 1 & 3) with metallicities Z = 0, 
10~ 6 and 1O _5 Z to compute the internal energy density 



and thermal pressure of the gas at a range of different 
densities. These values are used to construct look-up ta- 
bles which are then used by the SPH code to compute 
the pressure and internal energy at any required den- 
sity, for a given gas metallicity, via linear interpolation 
(in log-log space) between the tabulated values. The re- 
sulting equations of state are plotted in Figure Q] By 
using a tabulated equation of state, we avoid the large 
computational expense involved in solving the full ther- 
mal energy equation, while still obtaining qualitatively 
correct behavior. 

To model the star formation in this study we use 
sink particles, as described by Bate et al. (1995). This 
involves replacing the innermost parts of dense, self- 
gravitating regions of gas with particles that can both 
accrete further material from their surroundings and in- 
teract with other particles in the simulation via grav- 
ity. Sink particles are formed once an SPH particle and 
its neighbors are gravitationally bound, collapsing (nega- 
tive velocity divergence), and within an accretion radius, 
h acc , which is taken here to be 0.4 AU. Gravitational in- 
teractions between the sink particles and all other parti- 
cles in the simulation are also smoothed, self-consistently, 
to h acc . Our set-up allows us to identify sink particles 
as the direct progenitors of individual stars (for a mo re 
detailed discussion, see e.g. IWuchterl fc K lcsscn 2001). 

2.2. Setup and Initial conditions 

Our calculations are designed to start from the point 
where previous fragmentation calculations ended. It is 
now well established that the gas which falls into col- 
lapsing dark matter mini-halos undergoes a phase of 
fragmentation, resulting in the formation of large, self- 
gravitating clumps, with masses of order 100 M . The 
origin of this fragmentation is the rapid cooling driven by 
H2 that occurs at densities of around 10 - 10 4 cm -3 , and 
the masses of the clumps are similar to the Jeans mass 
at the corresponding density and temperature in this 
regime. We thus pick up the evolution fr om conditions 
simila r to those reached in the study of IBromm et al.l 
(|200ll) . Since the process of chemical enrichment involves 
supernovae events from the first stars we expect the gas 
to have a certain initial level of turbulence, which is nor- 
mally assumed to be absent during the formation of the 
very first stars from purely primordial gas. We also adopt 
a net angular momentum for the gas consistent with the 
results of simulations of cosmological structure forma- 
tion. 

The clouds in this study have a mass of 500 M , and 
are modelled using either 2 million or 25 million SPH 
particles. In the higher resolution calculations, this gives 
a particle mass of 2 x 10 -5 M and a mass resolution 
of 0.002 M (Bate & Burkert 1997), ro ughly 10 times 
smalle r than the opacity limit set by the lOmukai et al.l 
(2005) equation of state. These simulations therefore 
have roughly an order of magnitude of surplus resolu- 
tion. The low-resolution calculations have a mass resolu- 
tion roughly equal to the mass at which the gas becomes 
optically thick. To ensure that the Jeans condition is not 
violated, the sinks are always formed just before the opti- 
cally thick regime in all the simulations. Our clouds have 
an initial radius of 0.1 7pc, at an initial uniform density of 
5 x 10 5 cm~ 3 . This corresponds to an initial free-fall time 
of i ff = (37T/32G7T/9) 1 / 2 = 5.1 x 10 4 years. At this scale 
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Fig. 2. — Time evolution of the density distribution in the innermost 400 AU of the protogalactic halo shortly before and shortly after the 
formation of the first protostar at tgp. We plot only gas at densities above 10 10 cm -3 . The dynamical timescale at a density n = 10 13 cm -3 
is of the order of only 10 years. Dark dots indicate the location of protostars as identified by sink particles forming at n > 10 17 cm -3 . Note 
that without usage of sink particles we would not have been able to follow the build-up of the protostellar cluster beyond the formation of 
the first object. There are 177 protostars when we stop the calculation at t = r-gp + 420 yr. They occupy a region roughly a hundredth of 
the size of the initial cloud. With 18.7 Mq accreted at this stage, the stellar density is 2.25 X 10 9 Mg pc — 3 . 




Fig. 3. — We illustrate the onset of the fragmentation process in the high resolution Z = 10 — 5 Zq simulation. The graphs show the 
densities of the particles, plotted as a function of their x-position. Note that for each plot, the particle data has been centered on the region 
of interest. We show here results at three different output times, ranging from the time that the first star forms (t sf ) to 221 years afterwards. 
The densities lying between the two horizontal dashed lines denote the range over which dust cooling lowers the gas temperature. 



and density regime, the contributions from dark matter 
to the gravitational potential are small and are thus not 
taken into account in our computational set up. One 
can see from Figure Q] that the different gas metallicities 
have slightly different internal energies at the starting 
density. Thus, the Z < 10~ 6 Z Q calculations have an 
initial ratio of thermal to gravitational energy of a = 
-Etherm/I-Egravl =0.39, while the cooler Z = 10~ 5 Z Q cal- 



culations have an initial value of a = 0.32. All simu- 
lations are given a low level of initial turbulence, with 
the ratio of turbulent to gravitational potential energy 
-Eturb/I-Egravl = 0.1, and thus an RMS Mach number of 
M.RMS ~ 1. The clouds are set in initially uniform ro- 
tation, with (3 — E to t I 'l-Egrav] = 0.02. The conditions at 
the st art of our simulation are thus sim ilar to those at 
which iBromm. Coppi. fc Larson] (|2002f ) form their sink 
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Fig. 4. — Mass functions resulting from simulations with metallicities Z = 10 — 5 Zq (left-hand panel), Z = 10 -6 Zq (center panel), and 
Z = (right-hand panel). The plots refer to the point in each simulation at which 19 Mq of material has been accreted (which occurs at a 
slightly different time in each simulation). The mass resolutions are 0.002 Mq and 0.025 Mq for the high and low resolution simulations, 
respectively. Note the similarity between the results of the low-resolution and high-resolution simulations. The onset of dust-cooling in the 
Z = 10 — 5 Zq cloud results in a stellar cluster which has a mass function similar to that for present day stars, in that the majority of the 
mass resides in the lower-mass objects. This contrasts with the Z = 10 -6 Zq and primordial clouds, in which the bulk of the cluster mass 
is in high-mass stars. 



particles (see their Figure 5 for comparison) .We perform 
low-resolution simulations for all three metallicity cases, 
and high-resolution simulations for the primordial and 
Z = 1CT 5 Z Q cases. 

3. RESULTS 

We find that enrichment of the gas to a metallicity 
of only Z = 10~ 5 Z Q dramatically enhances fragmenta- 
tion within the dense, collapsing cloud. We first focus 
our discussion on the evolution of the Z = 10~ 5 Z Q gas 
cloud, before contrasting with the 10 -6 Zq and primor- 
dial clouds. 

The evolution of the high-resolution Z = 10 -5 Zq sim- 
ulation is illustrated in Figure [5J in a series of snapshots 
showing the density distribution of the gas. We show 
several stages in the collapse process, spanning a time 
interval from shortly before the formation of the first pro- 
tostar (as identified by the formation of a sink particle 
in the simulation) to 420 years afterwards. Note that we 
show only the innermost 1% of the full computational do- 
main. A complementary view is given in Figure El which 
shows the particle densities as a function of position. 

During the initial contraction, the cloud builds up a 
central core with a density of about n = 10 10 cm -3 . This 
core is supported by a combination of thermal pressure 
and rotation. Eventually, the core reaches high enough 
densities to go into free-fall collapse, and forms a single 
protostar. As more high angular momentum material 
falls to the center, the core evolves into a disk-like struc- 
ture with density inhomogeneities caused by low levels 
of turbulence. As it grows in mass, its density increases. 
When dust-induced cooling sets in, it fragments heavily 
into a tightly packed protostellar cluster within only a 
few hundred years. One can see this behavior in particle 
density-position plots in Figure [3l We stop the simu- 
lation 420 years after the formation of the first stellar 
object (sink particle). At this point, the core has formed 
177 stars. The evolution in the low-resolution simulation 
is very similar. The time between the formation of the 
first and second protostars is roughly 23 years, which is 
two orders of magnitude higher than the free-fall time at 



the density where the sinks are formed. Without the in- 
clusion of sink particles, we would only have been able to 
capture the formation of the first collapsing object which 
forms the first protostar: the formation of the accompa- 
nying cluster would have been missed entirely. 

The mass functions of the protostars at the end of the 
Z = 10~ 5 Zq simulations (both high and low resolution 
cases) are shown in Figure [4] (left-hand panel) . When 
our simulation is terminated, the sink particles hold ~ 
19 M of gas in total. The mass function peaks some- 
where below 0.1 M and ranges from below 0.01 M to 
about 5 Mq. It is important to stress here that this is 
not the final protostellar mass function. The continu- 
ing accretion of gas by the cluster will alter the mass 
function, as will mergers between the newly-formed pro- 
tostars (which cannot be followed using our current sink 
particle implementation). Protostellar feedback in the 
form of winds, jets and Hil regions may also play a role 
in determining the shape of the final stellar mass func- 
tion. However, a key point to note is that the chaotic 
evolution of a bound system such as this cluster ensures 
that a wide spread of stellar masses will persist. Some 
stars will enjoy favourable accretion at the expense of 
others that will be thrown out of the system (as can be 
seen in Figure \2§ , thus having their accretion effectively 
terminated (see for example, the discussions in Bonncll 
& Bate 2006 and Bonnell, Larson & Zinnecker 2007). 
The survival of some of the low mass stars formed in the 
cluster is therefore inevitable. 

Our ca lculations demonstra te that the dust-cooling 
model of lOmukai et al.l (|2005l ) can indeed lead to the 
formation of low-mass objects from gas with very low 
metallicity. This suggests that the transition from high- 
mass primordial stars to Population II stars with a more 
"normal" mass spectrum occurs early in the universe, 
at metallicities at or below Z w 10~ 5 Z Q . Our simula- 
tions even predict the existence of brown dwarfs in this 
metallicity range. Their discovery would be a critical 
test for our model of the formation of the first star clus- 
ter. The first hints of its validity come from the very 
low metallicity sub-giant stars that have recently been 
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Fig. 5. — The mass (upper panels) and the number of Jeans masses (lower panels) are both shown here as a function of gas density, 
with the high-resolution Z = 10 -5 Zq simulation in the left-hand plot and both the high-resolution primordial and low-resolution Z = 
10 -6 Zq simulations in the right-hand plot. In the upper panels, we plot the mass residing above a density, n gas , as a function of that 
density, as well as the Jeans mass. In the lower panels, we plot the number of Jeans masses as a function of density. In the left-hand 
panel, we show via the shaded areas, the heating (pink) and cooling (light blue) phases in the Z = 10 — 5 Zq EOS. In the right-hand plot, 
we show the conditions in the primordial (high-resolution) simulation (solid lines), and Z = 10 — 6 Zq simulation (dashed lines). The gas 
conditions are shown at the point of star formation and 50 years earlier, as well as at two instants after star formation, when the clouds 
have converted around 9.5 and 19 Mq of gas into protostars. Note that we only label the evolutionary stages in the lower panels, but 
the same progression with time, going from lower left-most lines to those at the top, applies also to the upper panels. 



discovered in the Galact ic halo (jChristlieb et al.1 [20021 : 
iBeers fc Christliebl 120051) . which have iron abundances 
less than 10 -5 times the solar value and masses below 
one solar mass, consistent with the range reported here. 

Turning our attention to the Z = and Z = 
10~ 6 Z Q calculations, we also find that fragmentation of 
the gas occurs, albeit at a much lower level than in the 
Z = 10~ 5 Z Q run. The mass functions from these sim- 
ulations are shown in Figure U (middle and right-hand 
panels), and are again taken when ~ 19 M of gas has 
been accreted onto the sink particles, the same amount 
as is accreted by the end of the Z = 10~ 5 Zq calculations. 

The primordial gas clouds form fewer protostars than 
the Z = 10~ 5 Zq clouds, with the high resolution sim- 
ulation forming 25 sink particles and the low-resolution 
simulation forming 22 sink particles. The mass functions 
are considerably flatter than the present day IMF, in 
agreement with the suggestion that Population III stars 
are typically very massive. The fragmentation in the 
Z = 10~ 6 Zq simulation is slightly more efficient than in 
the primordial case, with 33 objects forming. Again we 
stress that there is a delay of several local free-fall times 
between the formation of first and second protostars in 
these simulations: without the inclusion of sink particles, 
we would have missed the formation of the lower mass 
objects. 

4. CONDITIONS FOR FRAGMENTATION AND CLUSTER 
FORMATION 

We now discuss the physical origin of the fragmenta- 
tion in the Z = 10~ 5 Zq simulation, and investigate the 
properties of the forming cluster. First, we focus on the 
distribution of gas in the center of the halo right at the 



onset of star formation, i.e. at the time when we identify 
the first sink particle. Figure [5] shows (a) the distribution 
of rotational speed L/r relative to the Keplerian velocity 
Vkcp — (GM /r) 1 / 2 and (b) the specific angular momen- 
tum L in spherical mass shells around the halo center as 
a function of the enclosed mass. The blue curves denote 
the 10~ 5 Z Q case, while the red-dashed curves illustrate 
the behavior of zero metallicity g clS, clS discussed further 
below. It is evident that at the time when the first sink 
particle forms, rotational support plays only a minor role. 
The rotational velocity of gas in the inner parts of the 
halo is still sub-Keplerian by a factor or 4 to 5. This will 
change rapidly in the subsequent evolution as more and 
more higher angular momentum gas falls to the center. 
In (c) we plot the enclosed mass and in (d) the spherically 
averaged density as a function of distance from the cen- 
ter. The density roughly follows a power-law distribution 
with slope —2 ±0.1. 

We point out that the properties of the star forming 
clump in our model are virtually identical to those arising 
from full cosmological calculations taking into account 
the combined evolution of bary ons and dark matter over 
time (for comparison, see , e.g. lAbel. Bryan, fc Normanl 
120021 lYoshida et al\ 120061 ). If anything, one could ar- 
gue that the specific angular momentum we consider is 
somewhat on the low side. We therefore expect that 
fragmentation will also occur in full cosmological simula- 
tions, if these are continued bey ond the formation of the 
very first stellar object (see e.g. lYoshida et al.|[2007l for 
some preliminary evidence of this). 

To illustrate how the conditions for cluster formation 
arise we show in Figure[3]the density distribution perpen- 
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Fig. 6. — (a) Radially-averaged rotational velocity L/r normal- 
ized to the Keplerian velocity v^ cp = {GM cric (r) /r) 1 / 2 and (6) 
absolute value of the specific angular momentum L in spherical 
shells centered on the density peak as function of enclosed mass; 
(c) M onc (r) and (d) average local gas density n gas (r) as function 
of distance r to the center. Blue (solid line) curves denote gas with 
Z= 10 — 5 Zq and red (dashed line) curves denote gas with Z= 0. 
To guide your eye, a power-law slope of —2 is indicated in (d). 



dicular to the rotational axis of the system at different 
times. In Figure [5J in the top panel, we plot both the 
mass of gas which resides above a density n as well as the 
Jeans mass. In the bottom panel, we plot the number of 
Jeans masses, which is simply the mass at this density 
(or higher) divided by the corresponding Jean mass. 

The fragmentation in the low-metallicity gas (the Z = 
10~ 5 Z Q case) is the result of two key features in its ther- 
mal evolution. First, the rise in the EOS curve between 
densities 10 9 cm~ 3 and 10 11 cm -3 causes material to loiter 
at this point in the gravitational contraction. A similar 
behavior at densities around n — 10 3 cm -3 is discussed 
by Bromm et al. (2001). The rotationally stabilized disk- 
like structure, as seen in the plateau at n w 10 10 cm -3 
in Figure [31 is able to accumulate a significant amount 
of mass in this phase and only slowly increases in den- 
sity. Second, once the density exceeds n « 10 12 cm~ 3 , 
the sudden drop in the EOS curve lowers the critical 
mass for gravitational collapse by two orders of magni- 
tude. The Jeans mass in the gas at this stage is only 
Mj = 0.01 M Q , as visible in the top panel of Figure El 
The disk-like structure suddenly becomes highly unstable 
against gravitational collapse. We see this when looking 
at the behavior of the number of Jeans masses Nj above 
a certain density n in the bottom panel of FigureEJ which 
shortly after the formation of the first stellar object in- 
creases to values as high as JVj a 10 3 . Consequently, the 
disk fragments vigorously on timescales of several hun- 
dred years. A very dense cluster of embedded low-mass 
protostars builds up, and the protostars grow in mass by 
accretion from the available gas reservoir. The number 



of protostars formed by the end of our simulation (177) 
is nearly two orders of magnitude larger than the initial 
number of Jeans masses in the cloud set-up. 

Because the evolutionary timescale of the system is 
extremely short - the free-fall time at a density of n = 
10 13 cm -3 is of the order of 10 years - none of the pro- 
tostars that have formed by the time that we stop the 
simulation have yet commenced hydrogen burning, jus- 
tifying our decision to neglect the effects of protostellar 
feedback in this study. Heating of the dust due to the sig- 
nificant accretion luminosities of t he newly-formed pro- 
tostars will occur (|Krumholj|2005 l. but is unlikely to be 
important, as the temperature of the dust at the onset 
of dust-induced cooling is much higher than in a typical 
Galactic protostellar core (Td us t ~ 100 K or more, com- 
pared to ~ 10 K in the Galactic case). The rapid collapse 
and fragmentation of the gas also leav es no time for dy- 
namo amplification of magnetic fields (|Tan fe Blackmail! 
l200l . which in any case are expected to be weak and dy- 
namically unimportant i n primordial and very low metal- 
licity gas ()Widrowll2002l ). 

The cluster forming in our Z = 10~ 5 Zq simulation 
represents a very extreme analogue of the clustered star 
formation that we know dom inates in the present-day 
Universe (|Lada fc Ladall2003t) . A mere 420 years after 
the formation of the first object, the cluster has formed 
177 stars (see Figure [2]). These occupy a region of only 
around 400 AU, or 2 x 10~ 3 pc, in size, roughly a hun- 
dredth of the size of the initial cloud. With ~ 19 M Q 
accreted at this stage, the stellar density is 2.25 x 10 9 
M© pc -3 . This is about five orders of magnitude greater 
than the stellar density in the Tr apezium cluster in Orion 
(Hillenbra nd fc Hartmannlll998| ) and about a thousand 
times greater than that i n the core of 30 Doradus in the 
Large Magellanic Cloud (|Massev fc Hunter|[l99cl . This 
means that dynamical encounters will be extremely im- 
portant during the formation of the first star cluster. 
The violent environment causes stars to be thrown out 
of the denser regions of the cluster, slowing down their 
accretion. The stellar mass spectrum thus depends on 
both the details o f the initial fragment a tion process (e.g. 
as di scussed by Uappsen et al.l [2005t IClark fc Bonnell 
120051) as well as dynamical effects in the growing cluster 
(jBonnell et all 120011 : [Bonnell. Bate fc Vine Il2004h . This 
is different to present-day star formation, where the sit- 
uation is less clear-cut and the relative importance of 
these two processes may vary strongly f r om region to re- 
gion jKrumholz. McKee. fc Kleinl [20051: IBonnell fc Batd 



I2006t iBonnell. Larson fc Zinnecker ll2007l T 

Dynamical encounters in the extremely dense proto- 
cluster will also influence the binary fraction of the stars 
that form. Wide binaries will be rapidly disrupted, 
and so any binary sys tems that surv ive will be tightly- 
bound close binaries (jKroupal fl998). Recent observa- 
tions suggest that extremely metal-poor low-mass stars 
have a higher binary fraction than that found for nor- 
mal metal-poor stars, and that the period distribution 
of these binary systems is also skewed towards tight, 
short-period binaries (Lucatello et al. 2005). Mass trans- 
fer from a close binary companion may also be able to 
explain the extremely high [C/Fe] ratios measured in 
the most metal-poor star s currently known (jRvan et al.l 
120051 : iKomiva et aT1l2007t ). Our results suggest that these 
stars may originate in conditions similar to those that 
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Fig. 7. — Particle densities as a function of position in the low-resolution simulations, for the primordial (left), Z = 10 — 6 Zq (middle) 
and Z = 10~ 5 Zq simulations (right). The particles are plotted once the protostars in each simulation have accreted 19 Mq of gas. 



we find in our simulations. Further work aimed at im- 
proving our understanding of the binary statistics of low- 
metallicity stars formed by dust-induced fragmentation 
is clearly required. 

As mentioned in Section [3l the primordial and the 
Z = 10 -6 Zq cases als o exhibit som e fragm entation. 
Careful analysis of the lOmukai et~aTI (|2005[ ) EOS for 
zero-metallicity gas shows roughly isothermal behavior 
in the density range 10 14 cm~ 3 < n < 10 16 cm -3 , i.e. 
just before the gas becomes optically thick and begins 
to heat up adiabatically. Conservation of angular mo- 
mentum again leads to the build-up of a rotationally 
supported massive disk-like structure which then frag- 
ments into several objects. This is understandable, as 
isother mal disks are susceptible to gravitational insta- 
bility (|Bodenheimerl H995I ) once they have accumulated 
sufficient mass. Further, Goodwin et al. (2004a, b) show 
how even very low levels of turbulence can induce frag- 
mentation. Since turbulence creates local anisotropics in 
the angular momentum on all scales, it can always pro- 
vide some centrifugal support against gravitational col- 
lapse. This support can then provide a window in which 
fragmentation can occur. 

In addition to the quasi-isothermal behavior of the 
gas, both of the Z < 10~ 6 Z Q equations of state from 
Omukai ct al. (2005) also contain a brief phase of cool- 
ing during the collapse, which further aids fragmenta- 
tion. In the primordial case, this occurs very late in the 
collapse, just above n = 10 14 cm -3 , and is due to the on- 
set of efficient cooling f r om H2 collision- induced emission 
(|Omukai fc Nishill998t iRipamonti k. Abelll2004D . In the 
Z = 10~ 6 Z Q EOS, the cooling occurs earlier, at around 
n = 10 10 cm -3 , and is due to a combination of effects - 
enhanced H2 cooling resulting from the rapid increase in 
the molecular fraction at these densities due to efficient 
three-body H2 formation, and rotational and vibrational 
line cooling from H2O - that are present in the model of 
Omukai et al. (2005). One can see the emergence of struc- 
ture at these densities in the particle plots shown in Fig- 
ure [7] However, the effect is much less pronounced in the 
primordial case. Indeed a further low-resolution simula- 
tion of the primordial EOS in which this dip was removed 
yielded almost identical results, forming 17 sink particles 



instead of 22, suggesting that the quasi-isothermal na- 
ture of the gas is more important than this brief cooling 
phase. 

In comparison to the Z = 10~ 5 Z Q case, the strength 
of fragmentation in the Z < 10 -6 Z© calculations is 
weak, and only a few objects form for the combina- 
tion of total mass (M = 500 M Q ), angular momentum 
(E ro t/\E glav \ — 0.02), and level of initial turbulence 
(-Sturb/I-E'gravl = 0.1) that we adopted in our simula- 
tions. Consequently the stars (sink particles) forming in 
the Z < 10~ 6 Zq simulations are typically of higher mass 
than those in the Z = 10~ 5 Z © simulations, consistent 
with t he p redictions made by lAbel. Bryan, fc Normanl 
(12002ft and lBromm. Coppi. fc Larson! (|2002D . Although 
the recent high-re solution SPH simulat ion of primordial 
gas performed bv lYoshida et all (|2006f ), which predict a 
similar EOS to that used in our primordial simulations, 
find no fragmentation, they do not follow the evolution of 
the gas beyond the formation of the first collapsing core, 
since they do not include sink particles in their study, 
and so they may miss this subsequ ent phase of fragmen - 
tation. However, more recent work lYoshida et all (|2007l) . 
does show evidence of fragmentation on scales of around 
O.lpc. 

5. CAVEATS 

Although the equations of state taken from Omukai 
et al. (2005) provide an opportunity for vigorous frag- 
mentation at low metallicities, there are a number of 
questions which still need to be addressed. One imme- 
diate uncertainty is the applicability of the results de- 
rived from a 1-zone model to a full three-dimensional 
collapse calculation, which contains turbulence and thus 
local anisotropies in velocity and density structure. In 
particular, the existence of a strict relationship between 
temperature and density is unlikely, even when the cool- 
ing time of the gas is short compared to the dynamical 
time (Whitehouse & Bate 2006). If the gas were to col- 
lapse more slowly than is assumed in the Omukai et al. 
(2005) calculations, or if the opacity of the gas were to 
be lower, then less heating would occur prior to the onset 
of efficient dust cooling. 

On the other hand, if the collapse is faster than Omukai 
et al. (2005) assume or if the gas opacity is greater, then 



more heating would occur. Since the amount of heating 
that occurs during the loitering phase helps to accumu- 
late material at high densities, which then acts as a reser- 
voir for fragmentation, the amount of fragmentation will 
be sensitive to the thermal evolution of the gas during 
this phase of the collapse, and is therefore somewhat un- 
certain. However, we note that unless the temperature 
dip is eliminated entirely (which seems unlikely), frag- 
mentation will still occur along the lines outlined in this 
paper. We therefore believe our results to be qualita- 
tively, if not quantitatively, correct. 

A related issue is one of dust opacity. In their 
study, Omukai et al. (2005) use a dust model based 
on Pollack et al. (1994), which was designed for the 
study of Galactic dust. However, it is not at all clear 
that this is the appropriate model to use, given that 
high-redshift dust is likely to differ signific antly from 
local dust (see e.g. iTodini fc Ferraral 1200 ll ). Schnei- 
der et al. (2006) have performed a similar study to 
that of Omukai et al. ( 2005), but use a dust m odel 
based on the results of |Todini fc Ferraral (|2001f ) and 
ISchneider. Ferrara. fe Salvaterral (|2004D in place of the 
Pollack et al. (1994) dust model. They find that this 
leads to qualitatively similar behavior to the Omukai 
et al. (2005) model, but that the onset of effective dust 
cooling happens at lower metallicity, Z = 10 _6 Zq rather 
than Z = 10 -5 Z Q . However, the use of t his alternative 
dust m odel can also be questioned, since iNozawa et al.l 
(2003) have shown that the composition and proper- 
ties of the dust produced by population III supernovae 
are strongly dependent on the (poorly-constrained) de- 
gree of mixing assumed to occur in the ejecta, and since 
the model does not take into account the destruction of 
grains by thermal sp uttering in the reverse shoc k in the 
supernova remnant (|Bianchi fc S chneider! l200"7t ). or the 
growth of grains through accretio n or coagulation dur ing 
the protostellar collapse (see e.g. lFlower et al.~ll2005l ). 

Lastly, the total amount of heating produced by three 
body H2 formation is uncertain. This reaction, which 
takes place primarily at gas number densities between 
10 s and 10 13 cm -3 , is responsible for much of the rise in 
gas temperature prior to the dust-cooling phase. How- 
ever, there is a difference of two orders of magnitude be- 
tween the three-body H2 formation rate used in the sim- 
ulations of primordial star-formation performed by Abel 
et al. (2002) and the most recent theoretical determi- 
nation by Flower & Harris (2007), with other suggested 
rates spanning the full range in between. This introduces 
an additional uncertainty into the thermal evolution of 
the gas during the high-density loitering phase. 

6. CONCLUSIONS 

We have performed numerical simulations of star for- 
mation in very high density gas (in the range 5 x 10 5 < 
n < 10 16 cm" 3 ) in the early universe. The aim of the 
study was to investigate whether the dust-induced frag- 
mentation predicted by Omukai et al. (2005) and Schnei- 
der et al. (2002, 2006) does actually occur in realistic sys- 
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terns, and to begin to constrain the resulting IMF. The 
major differences of our work from the only previous nu- 
merical study of dust-induced fragmentation (Tsuribe & 
Omukai 2006) are the inclusion of the effects of rotation, 
which prove to be of vital importance, and the use of sink 
particles to capture the formation of multiple protostellar 
objects. 

Based on the equations of state reported by Omukai 
et al. (2005), our results show that fragmentation of 
protogalactic gas at very high densities above rt gas = 
10 12 cm -3 is almost unavoidable, as long as the angu- 
lar momentum is non-negligible. In this case, rotation 
leads to the build-up of a massive disk-like structure 
which provides the background for smaller-scale density 
fluctuations to grow, some of which become gravitation- 
ally unstable and collapse to form stars. At metallici- 
ties above Z ~ 10~ 5 Z Q dust cooling becomes effective 
at densities n gas ~ 10 12 cm" 3 and leads to a sudden 
drop of temperature which in turn induces vigorous frag- 
mentation. A very dense cluster of low-mass protostars 
builds up, which we refer to as the first stellar cluster. 
The mass spectrum peaks below 1 M© , w hich is similar 
to the value in the solar neighborhood (|Kroupal 120021 : 
IChabrierl l2003h and is also comparable to the mass of 
the very low metallicity subgiant stars recently discov- 
ered in the halo of our M ilky Way (jChristlieb et al.ll2002t 
iBeers fc Christ lieb 2005). If the dust induced cooling 
model proposed by Omukai et al (2005) is accurate, then 
the high-density, low-metallicity, fragmentation we de- 
scribe here may be the dominant process which shapes 
the stellar mass function. 

We find that even purely primordial Z = gas with suf- 
ficient rotation may fragment at densities 10 14 cm" 3 < 
n < 10 16 cm" 3 . In this density range, zero metallic- 
ity gas is roughly isothermal and the disk-like struc- 
ture that forms due to angular momentum conserva- 
tion is marginally unstable. It fragments into several, 
quite massive objects, thus supporting the hypothesis 
that metal-free stars should have masses in excess of sev- 
eral tens of Mq . Similar behavior is found for gas with 
a metallicity ofZ = 10" 6 Z Q . 
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